Submitted to ApJ. 



Carbon isotope fractionation in protoplanetary disks 

Paul M. Woods^ and Karen Willacy 

Jet Propulsion Laboratory, California Institute of Technology, 
4800 Oak Grove Drive, Pasadena, California, 91109, USA. 

paul . woodsOmanchester .ac.uk 
ABSTRACT 

We investigate the gas-phase and grain-surface chemistry in the inner 30 AU 
of a typical protoplanetary disk using a new model which calculates the gas 
temperature by solving the gas heating and cooling balance and which has an 
improved treatment of the UV radiation field. We discuss inner-disk chemistry 
in general, obtaining excellent agreement with recent observations which have 
probed the material in the inner regions of protoplanetary disks. We also apply 
our model to study the isotopic fractionation of carbon. Results show that the 
fractionation ratio, ^^C/^^C, of the system varies with radius and height in the 
disk. Different behaviour is seen in the fractionation of different species. We 
compare our results with ^^C/^^C ratios in the Solar System comets, and find a 
stark contrast, indicative of reprocessing. 

Subject headings: astrochemistry — solar system: formation — planetary sys- 
tems: protoplanetary disks 



1. Introduction 



Studying the chemistry of the inner regions of protoplanetary disks (PPDs) is a step 
into the unknown. Until recently observations of molecules in PPDs could not probe the 
warmer, inner regions of the disk where planets form. Detections of CO, HCO"*", H2CO 



C2H, CS, SiO, HNC, CN and HCN flDutrev et al.lll997l : iKastner et al.lll997l : lOi et al.ll2003 
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ester, Oxford Road, Manchester M13 9PL, UK 



- 2 - 



Thi et al.l l2004l : ISemenov et al.l |2005| ) at millimetre wavelengths can only tell us about the 
physics and chemistry in the cold conditions at radii greater than ~50AU because of the 
limits of millimetre-wave sensitivity and spatial resolution. We have to go to the infrared 
to investigate the hotter material, and there comparativ e ly few molecules have been seen. 
Initially CO JNaiita et allboosl : ISrittain et ZII2OO3I : Isiake fc Booeertl l2004h and H2O 
(ICarr et al.ll2004j ) were dete cted. Recently, the Spitzer Space Telesco pe has added OH, C2H2, 
CO2 and HON to the tally (ICarr fc Najitall2008l : iLahuis et al.ll2006l ) and a few detections of 
very large polycyclic aromatic hydrocarbon (PAH) molecules have been made (IGeers et al. 
20061). The Kec k Interferometer has also brought us detections of hot HCN and C2H2 



(iGibb et al.l 120071 ). The regions of the disk close to the star are of immense interest because 
they are where terrestrial (and larger) planets begin to form, as well as comets. Molecules 
in the inner 30 or 40 AU could be the building blocks of the kinds of pre-biotic and biotic 
molecules we see in the Solar System today. 

In the absence of observational evidence, we are left with theoretical work. This is 
valuable in its own right for understanding the processes in these regions, and for predict- 
ing observations for future proposals and even future technologies, such as ALMAl^, which 
will be able to probe these hidden inner disks. P revious modelling attempts of these very 



inner regions of PPDs (R<10AU) have been few flMarkwick et al.l l2002l: Millar et al.l 12003 



Ilgner et al.l 12004 : lllgner fc NelsonI l2006l : lAgundez et al.l |2008| ) , but have shown that these 
inner reg ions of PPDs are rich in molecules, including some complex molecules, such as 
benzene JWoods fc WillacvlboOTh . 



In this paper we present chemical models of a protoplanetary T Tauri-stage disk and a 
Minimum Mass Solar Nebula (MMSN), paying particular attention to the fractionation of 
carbon-bearing species. The fractionation in a particular molecule is defined as the abun- 
dance of the ^^C-bearing variant over the ^^C-bearing variant!^. Fractionation in carbon 
gives us a method of labelling: identifying, for instance, different regions of the disk and also 
identifying the different evolutionary processes involved in the formation of the disk - each 
process leaves its isotopic signature. Thus initially we study a typical protoplanetary disk, 
the physical aspects of which are explained in Sect. [2l We include a treatment of gas heating 
and cooling, which is expanded upon fully in Appendix |X1 In Sect. [3] we give a description 
of the chemical network, and the various chemical and photo-fractionation mechanisms in- 
cluding isotope-exchange reactions. Finally, we present our findings on carbon fractionation 
in protoplanetary disks, and set our results in both a Solar System and interstellar context. 



^The Atacama Large Millimeter Array, due for completion in 2012 (Iwww . alma^. inf op 

^When there is more than one ^'^C per molecule, the fractionation ratio is taken to be [^^C]/[^'^C] 
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2. The disk model: Physics 



In order to accurately model a PPD one has to take into account many physical and 
chemical processes. One needs a thorough understanding of disk hydrodynamics and magne- 
tohydrodynamics, turbulence, accretion, stellar spectrum, thermal balance, radiation trans- 
port, chemical composition, dust composition, etc. Whilst there are efforts to understand 
disks comprehensively and to model them accordingly with the appropriate feedback and 
interactions, such models are c omputationally very expensive and are very limited in their 
extent (e.g.. iTurner et al.ll2007l ). Thus simplifications are often made, effectively decoupling 
different elements of the situation. 

A major simplification is the separation of chemistry and physics when solving the 
equations of hydrodynamics for the disk temperature and density structure. In order to 
make the problem tractable a static hydrodynamic disk model is used to provide densities 
and dust tempe ratures throughout the di sk. From this the UV photon distribution can be 



calculated (e.g., Ivan Zadelhoff et al.ll2003l ). assuming contributions to the UV field from the 



stellar source and the interstellar medium. Few models have been developed th at deter 



mine the structure of t he disk and the radiative transfer self-consistently (e.g., iNomura 



2OO2I : iMillar et al.l l2003r). The final para meter required is that of the gas temperature. 



Some models (e.g., iMarkwick et al.l |2002| ) have assumed that the gas and dust tempera- 



tures ar e identical but this signifi cantly underestimates the gas temperature in the disk 
surface (IKamp &: Dullemondll2004l ). We calculate the gas temperature in the disk from the 
density and dust temperature profiles using a heating-cooling balance technique. This is a 
simpler approach than a finding a self-consistent solution (see lAikawa fc Nomural |2006| . for 
example), but neglects the feedback interaction between changes in the gas temperature and 
gas density. Once all the physical parameters have been determined they can be used as 
inputs for the chemical model. 

In reality, protoplanetary disks are not static: there is the large-scale motion of material 
accreting onto the star. There are also convective motions, often parameterised into vertical 
mixing motions and radial mixing motions. On all scales there is turbulence, which is mostly 
likely driven by magnetic fields in the disk. Turbulenc e drives mixing in both radial and 
vertical directions a r id models of niixing in the inner (jllgner et al.l 120041 ) and outer disk 
(IWillacy et al.ll2006l : ISemenov et al.l l2006l ) have shown this to be important. Material will 
also be gradually accreted radially towards the central star. Here we ignore turbulent mixing 
but mimic the accretion fiow of material moving toward the star by moving a number of 
parcels of gas from the outer edge of the model at 35 AU inward, along lines of constant 
scaleheight (zh). These parcels are spaced vertically in O.lzh intervals above the midplane 
of the disk, and we assume axisymmetry about the midplane. Each parcel will fiow into 
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the star in a time which can be derived from the principle of mass conservation. The radial 
velocity is given by: 

M 



with R, the radius, M, the accretion rate and S, the surface density of the disk, 
accretion timescale, At, is approximately AR/vj-. Thus, 

2n^RAR 



At 



M 



The 



(2) 



E sca les with 1/R fairly closely for R>10AU (ID'Alessio et al.ll200ll : [Andrews fc Williams 



20071 ). and so Eq. [2] only has a small dependence on R. A parcel starting at 35 AU will pass 



into the star in a time of approximately 0.41 Myr. In our framework of seventy radial points 
with a spacing of 0.5 AU, each parcel will spend ftiGOOOyr at each radial gridpoint before it 
is moved inwards. T his simplistic method has been treated in more detail, for instance, in 



Aikawa et al.l (119991 . eq. (10)). 



2.1. Disk structure 



To determine the physical struc ture of the disk we u s e a m odel kindly supplied to 



us by Paola D'Alessio and ba sed on 



D'Alessio et al l fll999l . I2OOII ). which models a flared 



disk using the a-formulation (jShakura fc Syunyaevi Il973l ). The model has a central star 
with a temperature (T^) of 4000 K, a mass (M^) of O- TM^d, a luminosity (L^) of 0.91©, 
and a radius (R^,) of 2.5 Rq (c.f., iGuUbring et al.lll998l ). The mass accretion rate (M) is 
10~^ yr~^ . The viscosity pa rameter, q;=0.01. The dust in the disk has an ISM size 
distribution (jPraine fc Ledll984j ) and is assumed to be well-mixed with the gas. The dust 
disk extends from ~0.1 AU to 300 AU and has a mass of 0.032 Mq. The model provides the 
density and grain temperature. The gas temperature is calculated separately, and details 
are given in the following section. 

From this model we have calculated a gas pressure scaleheight [zh) for the disk based 
on the distance over which the density drops to the 0.01% level. The radial dependence of 
Zfi can be fit by a power law: 

Zh = 0.0337i?^^/". (3) 
Our chemical calculations extend from the mi dp lane up to 6 Zh- 



2.2. Gas and dust temperature 
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Fig. 1. — Gas density and dust temperature profiles from tlie models of Paola D'Alessio, 
re-gridded to our 35x16 AU grid, spaced at 0.5 x 0.02 AU. 
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Fig. 2. — Results of the gas temperature determination for the inner 30 AU of the disk model. 
The slight "wobble" in the transition from the cold disk interior to the hot surface layers is 
due to a change from bare to mantled grains, which affects the heating due to H2 formation 
(Sect. ED. 
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Table 1. Gas heating and cooling mechanisms included. 



Mechanism 



Reference(s) 



Notes 



Heating 



Photoelectric effect 

H2 collisional de-excitation 
H2 photodissociation 
H2 formation 

C ionisation 
Cosmic rays 

Stellar X-rays 



Kamp fc van Zadelhoff (2001) 
Bakes fc Tielens (1994) 
Kamp fc van Zadelhoff (2001) 
Kamp fc van Zadelhoff (2001) 
Kamp fc DuUemond (2004) 
Cazaux fc Tielens f 2002b . 2004) 
Tielens fc HoUenbach (1985) 
Kamp fc van Zadelhoff (2001) 

Gorti fc HoUenbach (20041 
Shang et al. f2002l 



For silicate grains 

For small graphite and PAH grains 



H2 formation efficiency 

For S <150gcm~^ 
fUmebavashi fc Nakano 1981') 
X-ray ionisation rate 
Secondary effect 



Cooling 



Ol fine structure lines 



Cl fine structure lines 



Cll line at 157.7yum 

line of Ol at 6300 A 

H2 rovibrational lines 
CO rotational lines 



CH rotational lines 
Lyman-a line 
Gas-grain collisions 



Tielens k HoUenbach (1985) 
Kamp fc van Zadelhoff (2001) 
Tielens fc HoUenbach (1985) 
Kamp fc van Zadelhoff (^001) 
Tielens fc HoUenbach (1985) 
Kamp fc van Zadelhoff f200l1 
Tielens fc HoUenbach (1985) 
Sternberg fc Dalgarno (1989) 
Le Bourlot et al. (1999) 
Tielens fc HoUenbach (1985) 
Kamp fc van Zadelhoff f200l1 
Kamp fc van Zadelhoff f20011 
Kamp fc van Zadelhoff f20011 
Kamp fc van Zadelhoff f2001) 



Molecular line data 



Molecular line data 



Molecular line data 



Electron impact excitation 



Molecular line data 
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Until fairly recently PPD models assumed that the gas portion of the disk is in ther- 



mal equilibrium with the dust throughout its entirety (e.g.. iMarkwick et al.ll2002l ) . This 
can often be an underestimate of the gas temperature in regions of the disk close to the 
central radiation source, and especially in the tenuous surface which is exposed to the in- 
terstellar radiation field (ISRF), where conditions are similar to a photon-dominated region 
(PDR). This assumption c an underestimate the gas temperature by an order of magnitude 
faamp fc Dullemondl[20oi ). 



In an effort to more accurately model the isotope chemistry in PPDs, for which gas 
temperature is an important factor (see Sect. 13. ip . we have calculated the gas temperature 
separately from the dust temperature, by balancing heating and cooling terms. Several au- 
thors have previously taken this approach (IKamp &: van Zadelhofj|200ll : iKamp fc DuUemond 
2OO4J : iGorti fc Hollenbachll2004l : iTielens fc Hollenbachlll985h and we follow the prescriptions 
given in their papers, adapting them to T Tauri disks where necessary (see Table [1] for 
details). 

We take into account seven mechanisms which affect the gas heating rate in the disk 
- the photoelectric effect (Ppe), collisional de-excitation of H2 (Pcdx), photodissociation 
of H2 (Pphd), formation of H2 (rporm), ionisation of C (Fcion), cosmic rays (Pc-ray), and 
stellar X-rays (Px-ray)- Gas-grain collisions can also act as a heating mechanism when 
the dust temperature exceeds the gas temperature. However, since this rarely happens 
in the disk model, gas-grain collisions act mostly as a cooling mechanism (Ag_g)- Other 
cooling mechanisms included are: the fine-structure lines of atomic oxygen at 63.2, 145.6 
and 44.0/xm (Aqi), the fine-structure lines of atomic carbon at 609.2, 229.9 and 369.0/im 
(Aci), the 157.7/im line of singly-ionised carbon (Acii), the metastable line of atomic 
oxygen at 6300 A (Aoesoo); rovibrational lines of H2 (Ana), 25 rotational lines of CO (Aco), 
rotational cooling of the CH radical (Ach) and Lyman-a cooling (ALya). Further details are 
given in Appendix |X1 

Thus the temperature of the gas is calculated by balancing the heating and cooling 
rates. We use the mod i fied B rent's method of root-finding to solve for the gas temperature, 
as found in Press et al. Jl992h . This method works well in regions where the gas temperature 
function is continuous. Results of the gas temperature calculation are shown in Fig. [2J 

To increase computational efficiency we only perform this calculation once at every grid 
point, rather than iterating to a solution. Thus the gas temperature is calculated first using 
chemical abundances from the previous gridpoint and the new temperature used to calculate 
the chemical reaction rates for the current grid point. This is acceptable since there is only 
a weak coupling between the gas temperature calculation and the chemistry: any changes in 
temperature in response to a small change in chemistry will be small, and vice versa. 
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2.3. Radiative transfer 



Protoplanetary disks are subject to radiation from the central stellar source and from 
the ISRF. In order to calculate the UV radiation field at any o ne point in our disk mod el, we 
make use of the direct /diffuse treatment of FUV photons from lRichling &: Yorkd (120001 ) . We 
include the effects of radiation from both the central protostar and from the external UV field, 
and also include the effects of scatteri ng of photons by dust grain s. Absorption and scat tering 
coefficients are based on the work of iDraine fc Led (119841 ) and iPreibisch et al.l (119931 ). who 
assume a mixture of silicates, amorphous carbons and dirty ice-coated silicate grains. The 
disk is divided up into square cells (A2;=Ar=0.02 AU) and the radiative transfer equation 
is solved separately for both sources of radi ation, and for the direct and diffuse (scattered) 
components of each UV radiation field (see iRichling fc Yorkd l2000l . for details). The total 
UV fiux at a given point is the sum of the contributions from the stellar and interstellar UV 
fields, that is, we employ a I+ID approach rather than the full 2D treatment calculated by 
Richling fc Yorkel J2000h . 



The radiative transfer code is not able to account for the difference between the shape of 
the interstellar UV field and that generated by the spectrum of the young star. Observations 
have shown that T Tauri stellar fields are different to the interstellar radiation field, and can 
be dominated by strong emission features (such as Lyman-a). This can significantly impact 



the chemistry of the surface layers of the disk since 



such as OH and CH4, and not others. iBergin et al. (120031) found that t his effect can accoun t 



yma n-tt can dissociate some molecules. 



for high CN/HCN ratios observed in some disks (jPutrey et al.lll997l : iKastner et al.lll997l ). 
The strength of the T Tauri stellar radiation field has been estimated from FUSE observations 
as a few hundred times Gq at a radius of 100 AU (IBergin et al. I l2003h . where Gq is the standard 
interstellar radiation field strength. Here we adopt a value of 500 Gq at a radius of 100 AU. 

We also calculate pho toionisation by X-rays assuming an X-ray luminosity of 10~^ L^, 
(IGorti fc HollenbachI |2004| ) or 9 x 1O~^L0. The calculation of the X-ray photoionisation 
rate is described in more detail in Sec. 13.21 and the contribution of X-ray heating to the gas 
temperature is discussed in Appendix lA.ll 



2.4. Self-shielding 



Self-shielding by the abundant molecules H2 and CO moderates the effect of photodis- 
sociation. We incorporate this mechanis m in to our mo d el by using the "shielding factor" 
approach of Ivan Dishoeck fc Blackl (119881 ) and lLee et al.l (119961 ). For CO we use the data o f 
van Dishoeck fc Black Jl988l . see Sect. 13.1.21 for more detail), and for H2, iLee et al.l (119961 ). 
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These papers deal with interstellar clouds, for which the turbulent line width can be larger 
than the thermal l ine w idths in PPDs. We use the re-normalisation technique proposed by 
Aikawa &: HerbstI (Il999l ) to account for this difference, where the column density of H2 used 
in the self-shielding calculation is reduced by the factor C5/3kms~^ (where cs is the sound 
speed of the gas, and 3kms~^ is the turbulent line width assumed in the cited self-shielding 
calculations, for the ISM). In the inner disk, w here the temperatu r e and the sound speed are 
highest, this has less effect than that found by lAikawa fc HerbstI (119991 ) in the outer disk. 



We follow lAikawa fc HerbstI ( 1l999l ) in assuming that H2 and CO are only dissociated by 
the ISRF and not by stellar UV because of the difficulty of solving the equation of radiative 
transfer simultaneously in two dimensions. The stellar radiation is only effective in the 
relatively low density surface regions (above z/R~0.3, or 3.5 2;/^ at R=30AU), and thus has 
little effect on the chemistry in the molecular regions we wish to study. 



3. The disk model: Chemistry 

3.1. Isotopes of carbon 

Our knowledge of the chemistry of ^^C-bearing species stems from 30 years ago, and 
has advanced little since. This makes it difficult to accurately study a system which in- 
cludes ^^C-bearing species because reaction rates are simply not known. There is some data 
available on differences in the zero-point energy between ^^C- and ^^C-bearing species, but 
not for all spe cies, and the Kinetic Isotope Effect (KIE) has been used in theoretical work 



by some (e.g.. lYoung]|2006l ) to account for the differences in molecular mass, but this effect 
is complicated by the actual mechanism of reaction, i.e., the KIE only applies when the 
bond being broken/created is a bond between ^^C and another atom. However, there are 
two mechanisms which are known to fractionate carbon isotopes - fractionation through 
chemical exchanges, and through photodissociation. 



3.1.1. Fractionation through chemical exchange reactions 



Smith fc Adamsl (Il980l ) performed laboratory studies of isotope exchange reactions, and 



revealed that the most important reaction for exchanging carbon isotopes is: 

13Q+ ^12 QQ _13 QQ ^12 C++AE, 



(4) 



as predicted by IWatson et al.l (|l976l ). AE is the zero-point energy difference between the 
reactants and products, and is taken to be 35 K for this reaction. This small energy difference 
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makes the forward reacti on more efficie i it at low temperatures. The rate of this reaction 
was initially measured by IWatson et ahl (Il976[) at 300 K, and at 80, 200, 300 and 510 K by 



Smith fc Adamj fll980[) 



recent calc ulation bylLohrl ( 



anger et al. 



(119841 ) then calculated rates down to 5 K. A more 



with both ISmith fc Adamd ( 



99 



^ prod uced rates from 10-1 000 K, in reasonable agreement 



1980) and iLanger et al.l (119841 ) . The differenc e in the forward 



and b ackward rates of Eq. |l]is small, but it is still discernible even at 300 K (ISmith fc Adams 
1980h . 



Smith fc Adamd (1l980l ) also measured the rate of another isotope-exchange reaction at 



the same temperatures as before: 



H^^CO+ V CO ^ H^''CO+ V CO+AE 



13 



13, 



12 



(5) 



with AE=9K (ILanger et al.lll984l ). This is only effective at very low temperatures, with 
virtu ally no difference in forw a rd and backw a rd rat es at 300 K and only a small difference at 
80 K (ISnaith &: Adamslll980l ). ILanger et al.l (119841 ) calculated rates over the range 5-300 K 
and iLohrl (119981 ) made rate calculations over 10-1 000 K. 



Other exchange reactions similar to Eqs. (jlj) and may occur. However, a chemical 
reaction will almost certainly dominate over an exchange reaction such as reactions (jlj) or 
dS]) except when it is energetically unfavourable to do so (due to complicated structural 
rearrangement, for instance). In practice, many molecules react chemically with C"*", for 
example, rather than und ergo carbon exchan ge. An exception to this is CS, for whic h 
AE ^26K in the reaction JWatson et al.lll976h . Another is CH3 JPakarno fc Blackl EgTGl ) . 
Lohrl (119981 ) suggests that the reaction between HOC^ and CO might be important for 
isotope exchange since it moves in the opposite way to the exchange involving HCO"*" (Eq. [S]), 
preferenti ally putting ^ ^C into CO. The difference in zero point energies is small, 2.5 K. 
Similarly, iLangen (119921 ) suggested the importance of the possible exchange between C^ and 
CN, which might compete with the photodestruction of CN in the upper layers of disks. 
The difference in zero point energy is 34 K. However, none of the rates of these exchange 
reactions have been measured or calculated to our knowledge, and we do not include them 
in our model. 

For convenience, we have fit both the forward and backward reaction rates for reactions 
dl]) and ([5]) given in the literature with a smooth function in the standard form. Thus for 
reaction (HI), 



^for 

k 



3.3 X 10"^° I 

V300K 

kfor exp(— 35 K/T) 



-0.448 



3 —1 

cm s 



3 —1 
cm s 



(6) 
(7) 
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and for reaction 



^for 
h 



2.6 X 10 



-10 



T 



-0.277 



cm s 



^300 

fcforexp(-9K/T) cmV (9) 

These fits are in excellent agreement with both th e experimenta l resul ts of Smith fc Adamsl 
(Il980l . agree to within 7%) and the calculations of iLanger et al.l (119841 . 3%) and iLohn (11998 , 
8% above 50 K). 



3.1.2. Fractionation through photodissociation 

Carbon monoxide is dissociated mostly via line absorption and thus as the column 
density of carbon monoxide towards the emitter increases, so the absorption line saturates 
as photons are absorbed by the intervening material. At some point, a carbon monoxide 
molecule becomes shielded from dissociation by the presence of other molecules along the 
line of sight - not only molecules of the same kind, but also isotopologues and atomic and 
molecular hydrogen, whose dissociation bands may overlap. The degree of self-shielding 
depends on this overlap, and thus the shielding of ^'^CO, say, will depend on the column 
density of ^^CO and H2, as well as that of ^'^CO. One must also take into account the 
attenuation provided by dust grains. 



T his complex situation has been simplified by a number of authors for different situa- 



tions. 



an DishoeckVAN DISHOECK fc Blackl (11988h treat ^^CO explicitly, and therefore we 



use their treatment in our model. They have tabulated "shielding factors" which modify the 
photodissociation rate according to the factors detailed above. These shielding factors de- 
pend primarily on the column densities of ^^CO and H2, with the other factors (including the 
isotope ratio) remaining fixed. Since we calculate the column densit y of ^^CO in our model 



(and this differs from the fixed value of A^(^^CO)/A^(^^CO)=45 used in Ivan Dishoeck fc Black 
19881 ). we assume iV(^^CO)=45xiV(^^CO) for the purpose of calculating the shielding factor 



for ^^CO only. In general, iV(^^CO)/A^(^^CO) is less than 45, and so we may be slightly 
overestimating the contribution made to the self-shielding of ^^CO by ^^CO. 

For given column densities of ^^CO and H2, ^'^CO is shielded less than ^^CO, and thus 
is preferentially photodissociated. So as one proceeds from a region of high extinction to a 
region of low extinction, one would see ^^CO being photodissociated where ^^CO is already 
self-shielded, and hence an increase in the fractionation ratio, ^^CO/^^CO, i n this region. 



This effect can be seen observationally at the edges of molecular clouds (e.g.. ISheffer et al. 



19921 ). 
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This situation is complicated furth er by some recen t prehminary work showing that 
self-shielding is temperature dependent ( iLyons et al.ll2007l ). Since the predissociation bands 
for H2 and CO are thermally broadened and there is the possibility of overlap with adjacent 
electronic states, additional CO vibrational bands have to be considered at high tempera- 
tures. This work will be particularly applicable to the temperatures found in PPDs, and we 
look forward to the results of these investigations. 



3.1.3. Interstellar and Solar System context 



The isotope ratio for carbon (^^C/^^C) in the Solar System is widel y accepted to be 
89 JAnders fc Grevesselll98^ : Iciavton fc Nittlei][20oi : iMeibom et al.lboOTh. although recent 
measurements of the solar photosphere have indicated a ratio of 80±1 (lAyres et al.ll2006l ). 
This i s greater than in the local interstellar medium (ISM), where the value is taken to 
be 77 (IWilson fc Roodlll994l). greater than the Orion Bar region (^^C/^^C~60; iKeene et al. 



Tai,.,....^P...ia.«: aud much greater tharr the Galactic Cerrtre (-^ ^7^ 



Milam et al.ll2005l : iLanger et al.lll984l ) . This galactic gr adient (iLanger fc Penziaslll990l ) is due 
to the higher star formation rate in the inner Galaxy (jTosilll982l ). where th e fraction of ^^C 



has b een enhanced by the ^^C-rich ejecta of evolved, intermediate-mass stars (llben fc Renzini 
19831 ) in the time since the formation of the Solar System. 



Wilson fc Rood! (119941 ) give a numerical evaluation of how the isotope ratio changes 



with galactocent ric distance. Since the Sun formed approximately 1.9 kpc closer to the 
Gala ctic Centre ( Wielen et al.l ^9961 ) than its present location (7.94 kpc, Eisenhauer et al. 
20031 ). presumably it would h ave formed in a re g ion w ith a lower ^^C/^^C ratio, viz., Ri67 



ignoring temporal evolution. IWielen fc Wilson! (jl997f ) discuss a method of incorporating 
the temporal evolution of the ISM, and derive a value of 62 for the region and time period 
in which the Solar System condensed. The large difference between this value and the 
present Solar System value of 89 indicates that the Solar System must have become either 
significantly enriched in ^^C, or there must have been a significant depletion of ^^C. Recent 
studies of iron isotopes in the Solar System have indica ted that the Sun most likely formed 
close to one or more mass ive stars (jHester et al.1 120041 ). which produce ^^C in the triple-a 
reaction in their interiors ( Timmes et al.l 1995 ). These stars, which go on to form Type II 
supernovae, may have contaminated the solar protoplanetary nebula with ^^C-rich material 
during its formation. Altern ative explanations for this difference in isotope ratios are X-ray 
flares (Feigelson et al.ll2002l ). cloud mergers, orbital diffusion or radial gas streaming (see 
Milam et al.l l2005l . for further details). It seems unlikely that the increase in the ^^C/^^C 
ratio above the interstellar value is due to the processing of material once the Solar System 



-14- 



had been established (ILecluse et al.l 119981 ) . 

In light of these factors, we choose an initial ^^C/^'^C ratio of 77, between the values 
of 62 and 89, in line with that of the present-day ISM ratio. The value of this ratio is not 
crucial to the chemistry, and similar fractionation levels are obtained for the entire range of 
values, 62-89. 



3.2. The reaction set 



The chemical reaction network is based on the UMIST99 gas-phase ratefile (iLe Teuff et al. 



20001 ). with two additions: 



1) we include grain surface reactions from lHasegawa fc HerbstI (119931 ) and lHasegawa et al. 

(I1992I ). We consider cosmi c ray h e ating and th e rmal desor ption, with update d binding en 



ergies from Bisschop et al. (2006), Oberg et al. ( 2005[) and 



Fraser et al. 



H, we use updated binding energies from ICazaux fc TielensI J2004j . l2002bh . who treat ph- 



(IgOOlh . For atomic 



ysisorption and chemisorption onto grains distinctly, with binding energies ii^Hp=600K and 
i?Hc=10 000K, respectively. In our calculations we assume that the sticking coefficient of H 
is 0.4, and that for all other species is 0.3. 



2) we include a simple treatment of X-ray chemistry similar to that by lGorti fc HoUenbach 

( I2OO4J ). i.e., we calculate an ionisation rate per atom, which depends on photon flux and 



cross-section: 



= 6.25 X 10^ 



10 



0.5 



(10) 



for a tom i. Fits to X-ray cross- sect ions (ux) for astrophysically relevant molecules are taken 
from Verner &: Yakovlev ( 1995 1^. F(E) is the X-ray photon flux at a radius R, as a function 
of energy 0.5< E <10keV. tx{E) = Nq(7x{E) is the X-ray extinction due to a column 
density A^o- We assume that X-ray ionisation leads to the loss of a single election, and 
that the ionisation rate for molecules is the sum of the rates for the constituent atoms. See 



Aikawa &: HerbstI (1200 ll ) for a more detailed treatment of X-ray chemistry. 



In total, our reaction network comprises 475 gas and grain species, and over 8 000 gas- 
phase and surface reactions. More than three-quarters of these rea ctions involve ^^C. There 
have been recent observations (ISakai et al.ll2007l : lTakano et al.lll998l ) of ^^C isotopologues and 
isotopomers which show that there may be differing properties for carbon atoms attached in 
different places in the molecule. However, due to a dearth of experimental data on how ^^C- 



■^Please contact the authors for our fitting coefficients 



- 15 - 



Table 2. Initial fractional abundances for abundant carbon-bearing species. These 
abundances are the result of 1 Myr in a molecular cloud model 



Species 


Gas phase 




Solid phase 


g-i^CX/g-i^CX 




1.7 (-7) 


169 


negligible 


- 


C3H2 


4.1 (-7) 


160 


1.9 (-6) 


133 


C2H2 


3.6 (-7) 


153 


1.2 (-6) 


141 


C4H 


5.3 (-7) 


152 


negligible 


- 


HCN 


2.4 (-8) 


152 


3.7 (-6) 


Ill 


HC3N 


1.1 (-7) 


150 


5.8 (-7) 


134 


CN 


1.3 (-7) 


146 


negligible 


- 


HNC 


1.0 (-8) 


144 


2.6 (-7) 


106 


C4 


7.2 (-8) 


t A A 

144 


1 * *1 1 

neghgible 




Ca 


9.8 (-8) 


134 


negligible 




C3H 


5.3 (-7) 


133 


negligible 




C4H2 


9.6 (-9) 


133 


8.5 (-7) 


136 


C2H 


6.5 (-8) 


118 


negligible 




CH 


2.4 (-8) 


113 


negligible 




CH3 


4.8 (-8) 


107 


negligible 




H2CO 


4.7 (-8) 


107 


2.1 (-7) 


106 


CH4 


2.3 (-6) 


102 


9.2 (-6) 


99 


C 


1.1 (-7) 


93 


negligible 




CO 


3.1 (-5) 


55 


2.6 (-6) 


55 


C02 


1.7 (-7) 


51 


6.4 (-7) 


55 



Note. — X {—y) here represents x x 10~^. Isotope ratios for species with more 
than one carbon atom are calculated by dividing the total abundance of ^^C by 
the total abundance of ^^C in that particular molecule. Other species of note, and 
their initial fractional abundances: H2 5.0x10"^, He 1.4xl0~^, g-H20 1.4x10"^, 
H 2.5xl0-^ g-NH3 1.6xl0-^ O 1.4xl0-^ OH 6.4x10"^. 
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bearing species react in comparison to ^^C-bearing species, we have assumed that reactions 
involving ^^C proceed at the same rate as their ^^C counterparts. We have taken account of 
the increased number of reaction products which comes from the inclusion of isotopic species. 
For example, the reaction: 

C2 + H — ^CH + C (11) 

with rate k now becomes: 

120^2(3 + H — > I'CH+i'C (12) 
^^C^^C + H — > ^^CH+^^C (13) 
12C13C + H (14) 
^^C^-'^C + H — y ^^CH+^^C, (15) 

with rates k, k/2, k/2 and k, respectively. We preserve functional groups such that, e.g., 

^^CH^2CN + He+ — > ^^CN CRj + Re (16) 
^3CH^2CN + He+ -h ^^CN CH+ + He (17) 

and we preserve double bonds in preference to single bonds: 

H^^CO+^^CH — > H^=^C0+^^CH2 (18) 
H^^^CO+^^cH ^ (19) 



As initial fractional abundances we use the outputs of an interstellar (IS) cloud model 
(r;,=2xl0^ cm~^, T=10K, 74v=10) which uses the same chemical network, and is allowed 
to run for 10^ yr. The inputs to this cloud model are the "low metal abundances "of 
Graedel et aDJl982h . viz. H:He:0:C:N:Si are 1:0.14:1.76x10-^:7.30x10-^:2.14x10-^:2.00x10- 
Our IS cloud mod e l is a simple single-point approximation, but it reproduces the results of 
Langer fc Graedell (119891 ) very well. We are in agreement to within a factor of 2 for impor- 
tant species such as CO, IICO+, O, CH, CII2 and H2CO, and a factor of five agreement with 

C2 and H2O. There is a lesser agreement with CN, HCN, C2ll, CH4 and OH due to the 
advances in the accuracy of reaction rate determinati on in the last 25 years. For instance, 
our model produces an overabundance (compared to iLanger fc Graedell Il989l ) of CN by a 
factor of ~20, due to a faster rate for the reaction: 



CH + N 



CN + H. 



(20) 



from 



Langer fc Graedell (119891) u se a rate of 4.5x10 



-12 

cm s 



Le Teuff et al.l (I2OOOI ) is 2.1x10 



-lOpj^Sg 1 



^ at 20 K, whereas the revised rate 
47 times faster. See Table [2| for input 
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abundances for select species, and also the fractionation ratios after the molecular cloud 
stage. 

The chemical frac tion ation which occurs i n an interstellar clo ud is explore d in d etail 
by iLanger et al.l (119841 ) and iLanger fc Graedell (119891 ) . In summary, iLanger et al.l (119841 ) are 
able to classify C-bearing species into three families - CO, HCO"'" and the "carbon isotope 
pool" - with distinct isotopic behaviours. The fractionation in CO is driven mainly by 
reaction (jl]), which favours the production of ^^CO at the low temperatures and densities 
of interstellar clouds, driving the ^^CO/^^CO ratio down. The fractionation of HCO"*" is 
driven both by reaction ([5]), which preferentially puts ^^C into HCO"*" at low temperatures, 
and also by its formation from elements of the carbon isotope pool, which are ^^C-enriched. 
The fractionation in these remaining carbon-bearing molecules is driven to high ^^C/^'^C 
ratios since the chemistry is based on C"*", for which ^^C^ is favoured at low temperatures 
(reaction H]) and much ^'^C is taken up in ^'^CO. 



Results 



4.1. Gas heating and cooling 



The gas temperature in our model of a circumstellar disk can reach ~8 000 K at the 
surface, vastly exceeding the temperature of the dust from the input dust models, which 
is a few hundred degrees in the same region. Such a n eff ect has also been fou nd by 
Kamp fc van Zadelhofj (|200ll ): lKamp fc DuUemondl kooj ) and loiassgold et al.l toQ^ using 
different disk models. Clearly there is a strong need for disk models which calculate dust 
and gas temperatures individually. The gas and dust are not well coupled coUisionally (and 
hence not in thermal equilibrium) in the disk surface, where the optical depth to the ISRF 
(r) is less than 0.5 (~3 2;/i). 

Heating by the photoelectric effect on small carbon and PAH grains is very important 
in the surface layers of the disk. This mechanism dominates almost all others in the upper 
third of the disk; gas heating due to H2 formation becomes the most effective mechanism 
in a limited region where H2 is photodissociated, around z=4-5zh- Intermediate layers of 
the disk are heated mainly by X-ray heating for radii >10AU, and by gas-grain collisions 
inside of this radius. The midplane is heated predominantly by coUisional de-excitation of 
H2 molecules. The lower two-thirds of the disk are cooled mostly by molecules - in the 
most part, CO molecules are the most effective coolant, although CH molecules cool the 
midplane very efficiently at larger radii, R>17AU. In the upper third of the disk, cooling 
by the forbidden lines of Ol is very effective, with Lyman-a cooling dominating at the very 
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surface of the disk, inside of 20 AU. 

The transition between the hot surface layers and the cool bulk of the disk shows some 
small oscillations in temperature (Fig. [2]). This is a result of the inclusion of icy grains, 
with a different H atom binding energy compared to bare grains in our model. On icy 
grains H atoms can only physisorb and therefore have a relatively short residence time. The 
oscillations occur in the region where the transition from bare to icy grains occur. 



4.2. Fractionation of carbon isotopes 

The disk can can be split into different regions depending on temperature, which is 
the dominating factor in the carbon fractionation of different species. The cold midplane 
region ("the cold region", T=30-100K) of the disk extends up to ^^3.5 Zh, which equates 
to 9AU at a radius of 30 AU. In this region the majority of molecules are frozen out onto 
the surface of dust grains; only very volatile species (CO, CH4) are in the gas phase. The 
warmer temperatures in the very inner part of the disk means that the ice line for water falls 
close to 2 AU. 

Above the cold region is a transition region ("the transition region", T=100-2 500 K), 
which generally has a thickness of ~0.8z/i (2AU at /2=30AU). The increase in temperature 
in this region causes molecules to evaporate from grain surfaces - H2O is one of the last species 
to persist on grains. UV extinction is low enough for most molecules to be photoprocessed, 
although H2, ^^CO and ^^CO remain shielded. 

Lastly, the surface region (above ~4.3 Zh) of the disk is heavily ionised and very hot 
("the hot region", T=2 500-8 000 K). Ionised species have a fractional abundance of ~10~^ 
in the surface region, indicating that most hydrogen is atomic and not yet ionised. These 
three regions are evident in Figs. [3] and [51 We will look at some characteristic species in 
detail and their properties in these three regions. 



4.2.1. CO, HCO+ and CO2 



CO is an important species in disks since it is the dominant gas-phase molecule other 
than H2 and can be used observationally to trace the bulk gas of the disk. It is also im- 
portant because it can be u sed to trace the vert ical temperature struct ure in disks, as in 
Guilloteau fc DutrevI Jl998h : lOartois et all J2003h and lPietu et"aD J2007h . 



Of the CO introduced into the disk from the interstellar medium at 35 AU in our model, 
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92% of that is in the gas phase, an d 8% is in the form of CO frozen onto dust grains. CO 
desorbs at low temperatures (26 K; lOberg et aLll2005l ) and in the inner disk temperatures 
are always above this, so CO quickly desorbs and very little ('d%) CO ice remains. Hence 
there is a high abundance of gaseous CO available to take part in chemical reactions. 

The fractionation ratio of CO in the disk has a limited range, varying from 25-77 
throughout the disk. This contrasts with observed ratios in diffuse inte rstellar clouds where 
there is a much wider range, 15<N(^^CO)/N(^^CO)<170 flLisztl 120071 1. The fractionation 
ratio of CO in the midplane does not change a great deal, 47-55 between 1 and 30 AU 
(Fig. [3]). However its high abundance means it is able to be involved in exchange reactions 
thereby altering the fractionation of other molecules, e.g., HCO"*", which shows an increase in 
fractionation with decreasing radius. This is a consequence of exchange reaction Under 
the high density conditions of the disk, this process is roughly in chemical equilibrium: 



r2(Hi2CO+) 



exp(-9K/T) 



CO) 
^^CO) 



(21) 
(22) 



n(Hi3C0+) ' 'n(i3C0) 

using the relation in Eq. (Q. Thus at temperatures of 32 K, where the ^^CO/^'^CO ratio is 55, 
we should expect H^^C0"'"/H^^C0"'"=41, and at temperatures of 270 K, 45, as seen. Thus at 
higher temperatures (T>60K), the fractionation of HCO^ begins to trace the fractionation 
of CO, at least for the midplane region of the disk where CO is the only source of HCO"*". The 
fractionation of HCO"*" is shown in Fig. HI Within 10 AU, HCO"*" reacts with increasingly 
abundant hydrocarbons (e.g., C3H4, C4H2) and ammonia to create reactive ions and to 
reform CO, and within 8AU these reactions become faster than the interchange between 
CO and HCO+. 

At radii of a few AU, the reaction between CO and OH is the major contributor to the 
formation of CO2. This reaction has a moderate activation barrier, and thus proceeds more 
rapidly at higher temperatures. This reaction, which is important for both atmospheric 
and combustion chemistry, i s one of the few for A vhich mass-independent isotope effects 



have been investigated (c.f., IChen fc MarcusI l2005l : [Stephens et al. 



1980; Smit et al. 1982 



Roeckmann et al.lll998l ). These calculations and experiments find a small difference in the 
reaction rates of ^^CO and ^^CO with OH, but they were carried out at significantly higher 
pressures than are found in protoplanetary disks. At the lowest pressures considered, they 
found that the reaction involving ^'^CO is faster than that involving ^^CO by less than 1%. 
Given that this percentage is dwarfed by the uncertainties in published reaction rates, the 
Kinetic Isotope Effect (of which this is an example) is unlikely to affect our results. 

The fractionation ratio of CO ice in the midplane resembles that of gaseous CO. CO2 
ice, with a higher binding energy than CO ice, has a fractionation ratio which decreases from 
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Fig. 3. — Carbon monoxide fractionation shown as the ratio of ^^CO/^^CO throughout the 
disk. Three different regions can be seen in the fractionation data - the hot region at the top 
of the disk, the cold region at the bottom and the transition region in between. The sohd 
hues in this and following figures represent levels of optical depth to the ISRF, r=0.1, 1, 10. 
The dotted line indicates the designated surface of the disk at 6 Zh- 




5 10 15 20 25 30 5 10 15 20 25 30 

Radius [AU] Radius [AU] 



Fig. 4. — The fractionation and distribution of the formyl ion, HCO"'". White areas in this 
and subsequent figures indicate regions of negligible fractional abundance, a;(X)<10~^^. 
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its input value of 55 to a value slightly higher than that of CO (46) in the very inner part 
of the disk, 49. The value o f ^^CO^/^^COo ice ra tio has been determined to be 81±11 in 
a young protostar, Elias 29 (IBoogert et al.ll2000al ). This is unusually high, especially since 



the nascent cloud, p Op h, shows ^^C/^^C ratios typical of the local ISM (ICasassus et al. 



20051 : iBensch et al.ll200ll ). This value is also somewhat hi gher than the ^^CO/^^CO ice ratio 
along the same line of sight, 71±15 (IBoogert et al.ll2002l ). A suggestion for this high value 
( Boogert et al.ll2000bl ) is that CO2 may have been formed from C*^"'"'' rather than CO, as 
gene rally assumed. In our model, th e vast majority of CO2 ice comes from CO, in agreement 
with lEhrenfreund fc Schuttd (120001 ). Whatever the formation route, it seems that CO2 ice 
ratios are determined in the parent cloud rather than the hot core around the protostar 
dCharnlev eraDl2004h . 



Towards the top of the cold region in the disk {z ~2.7zft) there is a layer of low abun- 
dance of CO. Above this layer, CO can be destroyed by reactions with He"*", from the X-ray 
ionisation of He. Some carbon ions resulting from this destruction end up forming hydrocar- 
bons on the surfaces of dust grains. At this height in the disk, the grain temperature is just 
high enough that hydrocarbons can thermally desorb, ensuring that the carbon contained in 
them is not lost from the gas. In contrast, in the low abundance layer at 2.7 z^, hydrocarbons 
that form on the grains are retained, resulting in a loss of carbon from the gas, made evident 
in a loss of CO. Below this low abundance layer. X-rays do not penetrate and the destruction 
rate of CO by reactions with He"*" or B.^ is significantly lower. This effect is greater in the 
MMSN model (Sect. [H]) because that model has a higher density and a lower temperature, 
i.e., more efficient freezeout and less efficient desorption for a given column density. 

In the transition region of the disk, fractionation is driven strongly by photoprocesses. 
HCO"*" becomes abundant, with fractional abundances of 1-10x10""^°, and consequently 
increases the fractionation ratio of CO, through reaction (jS]). Due to self-shielding effects, CO 
is more resilient to photodissociation than other carbon-bearing species and thus is molecular 
in a region where other molecules are starting to become dissociated by UV photons. There 
is also a difference in the degree of self-shielding between ^^CO and ^^CO, and thus there is a 
narrow layer in which ^^CO is photodissociated and ^^CO is not. This "photo-fractionation" 
layer occurs at the very top of the transition region {z ~ 4.1 Zh), and is minimal in thickness 
due to the low column densities to the dissociating UV radiation fields at the top of the 
disk. This region of the disk is also the one in which the exchange between CO and C"*" 
(reaction H]) is most important. C^ and ^'^C^ are abundant because CO and ^^CO are being 
photodestroyed, and isotope exchange and photodissociation compete. 
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Fig. 5. — Variation of ^^C+/^^C+ and ^^C/^^C throughout the disk. 
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Fractionation in the top layer of the disk is relatively easy to understand, since all but 
trace amounts of carbon are to be found in the form of C"*" (and ^'^C^). Thus it follows that 
the fractionation ratio in this region must reflect the "input" value. For our fiducial model, 
this input value is 77, and Fig. O shows this value in the upper layer of the disk. 0+ is 
produced very rapidly, on a timescale of days, mainly by photoionisation of C. This is much 
faster than a vertical or radial mixing timescale, and thus the fractionation in the upper 
ionised region is likely to persist as a long-lived feature of the disk. This provides us with a 
simple way to quantify the carbon fractionation in an observed disk. 

The fractionation pattern in the upper, hot layer of the disk shows the need for gas 
temperature calculations in modelling. Assuming that gas and dust temperatures are iden- 
tical produces a very uniform disk, with slight increases in fractionation in the transition 
region, where the difference in self-shielding factors for carbon monoxide is evident. The 
upper region of the disk has the same degree of fractionation as lower levels, and in general, 
carbon-bearing species do not reflect the input fractionation ratio. In our model, where the 
gas temperature is calculated separately from the dust temperature, the activation barrier 
of reaction (jlj) becomes negligible at sufficiently high temperatures, such as those found in 
the disk surface layers. Thus the fractionations in and CO are averaged into the input 
fractionation ratio of 77. 



4-2.3. H2CO, C and the carbon isotope pool 
The fractionation i n H^CQ in interstel lar clouds provides one of the upper bounds to 



the total ^^C/^'^C ratio (ILanger et al.lll984j ). However, this is not the case in disks, where 
the fractionation in H2CO varies from ~60-100 (FiglH]). H2CO is formed almost exclusively 
by the reaction between CH3 and O in the gas phase, and thus the fractionation in H2CO 
reflects that in atomic carbon. 

The fractionation ratio of atomic carbon in the disk varies from 16-110 (see FigE]). In 
general this can be further constrained to ~45-110, with a thin layer 3.3-4.5 Zh (0.5<R<10 AU) 
where the ratio drops precipitously to ~15. This drop is due to a slight enhancement of the 
photodissociation rate of ^^CO over ^^CO due to the differences in self-shielding, causing 
the abundance of ^^C to rise. Atomic carbon is the basis for the formation of many hydro- 
carbons, and thus species in the carbon isotope pool (e.g., CH, CH4, etc.) will follow the 
fractionation in C. 
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4- 2.4- Nitrogen-bearing species 

Solid HCN is a major repository of carbon in the cold regions of the disk, storing up 
to 5% of all available carbon, as well as ~20% of all nitrogen. Desorption of HCN from 
grain surfaces becomes efficient at a height of ~3.2-3.4^/,. HCN is then destroyed in the gas 
phase by photons. What happens to the liberated nitrogen depends on position in the disk. 
At 3.5zh it is cycled back into HCN through reactions with carbon clusters (C„). At 3.7zh, 
the availability of oxygen-bearing reactants derived from thermally desorbed H2O is much 
greater than at 3.5zh, and so N is cycled into NO, and C is cycled into CO through reactions 
with OH. This causes a slight increase in the abundance of CO, which can be seen in the 
right panel of Fig. [31 

Solid HCN retains the fractionation ratio of the interstellar cloud model, H^^CN/H^^CN 
= 111. In the gas phase, in intermediate layers of the disk, the fractionation of HCN, HNC 
and CN is controlled by, and mimics, the fractionation of C. 



Discussion 



5.1. Comparison with previous disk models 



The most relevant work on inner disk chemistry is that by iMarkwick et al.l (120021 ) , who 
study the inner 10 AU of a static disk and assume that the gas temperature is equal to the 
dust temperature. Their physical model d iffers greatly from ours, although their chemical 
network is similar. iMarkwick et al.l (120021 ) adopt an accretion rate which is a factor of ten 
larger than ours - this will cause their surface density to be higher, and also their disk to be 
warmer, in general. The higher surface density means that fewer photons will penetrate the 
disk, and thus the ionisation profile of the disk is different (their Fig. 3 shows a significantly 
lower abundance of H+ - x(H"'")<10^^^ - than our model - x(H'*')<10~^). Furthermore, the 



Markwick et al.l (120021 ) model involves a 1 Mq star, 43% more massive than in our model. A 
higher stellar mass implies a greater stellar gravity and thus a thinner disk, which receives 
less stellar UV radiation. Their model also has a temperature inversion, which means that, 
surprisingly, they find that the majority of species are adsorbed onto dust grains in the 
surface of the disk at 10 AU. Given this, it is not surprising that the predictions of the two 
models are very different. Ionisation products such as HCO^ are much more abundant in our 
model (compare Fig. IHwith their Fig. 4), and volatile species are available in the gas phase 
at different heights and radii due to the difference in the dust temperature profile caused by 
the different accretion rate. 
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A more recent model from lAgundez et al.l (120081 ) treats the chemistry in the high den- 
sity PDR-hke regions of protoplanetary disks - the surface of the disk and the inner 3AU. 



Agundez et aLl (120081 ) use a solely gas-phase model to investigate the formation of small 



species such as HCN, C2H2 and CH4, introducing reactions with significant activation energy 
barriers (~1 400 K) which are not included in chemical networks based on int erstellar chem- 



i stry, such as the UMIST and Ohio State networks. The disk model used by lAgundez et al. 



( I2OO8I ) is similar to the one used here, and hence comparison is straightforward. Since they 
only consider the surface region of the disk, they calculate column densities vertically from 
a height Zi^, the height in the disk where the gas is well-shielded to UV radiation and where 
all the carbon is as CO. In our disk, this corresponds to a height of ~3 z^. Ag reement is very 



good for small species such as CO, CO2, CH4 and C2H2, although the model of lAgundez et al. 



(I2OO8I ) is deficient in H2O, OH, NH3 and HCN compared to our model. These differences may 
be explained by different photodissociation rates and activation barriers for these species. 
For instance, t h e acti vation barrier for the reaction H2 -|- CN — > HCN + H is 1 200 K in 



Agundez et al.l (120081 ). but only 820 K in this work, with this value taken from the NIST 
database. 



5.2. Significance for observations 



Observations of protoplanetary disks at long wavelengths (sub-millimetre and millime- 
tre) are sensitive to the cold gas of the midplane and outer regions (R>50-100 AU). Recent 
observations at infrared wavelengths are sensitive to warmer gas and dust, and thus have 
been able to probe within a few AU of the star. Future initiatives (e.g., ALMA) which 
give us high spatial resolution, should enable us to detect material at the radii of Earth-like 
planet-forming regions in other systems. 

Some recent observations have led to derivations of the ^^C/^'^C ratio in one edge-on (GV 
Tau) and one almost face-on (HL Tau) T Tauri disks, similar to the one we model. GV Tau is 
a binary system in whi ch the T Tauri star has a ste llar mass and effective temperature similar 
to that which we use (|White fc Hillenbrandl |2004| ) . The mass of the disk is estimated to be 



0.01 Mq ( Hogerheijde et alll998h. The ac"cretion rate is 20 times greater, 2x10 "^MqYi 



and the luminosity nearly 2 Lq (IWhite fc Hillenbrandll2004l ). HL Tau has a similar accretion 
rate and stellar temperature to GV Tau, but is nearly twice as massive and slightly less 
luminous (IWhite fc Hillenbrandl 120041). The di sk around HL Tau extends to at least 200 AU 
and has a mass of ^O.IM© (iGibb et al.ll2004l ). Thus there is a good basis for comparison 
between our model and these objects. 



Through a combination of fundamental and overtone lines of ^^CO and CO. IGibb et al. 
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Fig. 6. — Formaldehyde fractionation shown as the ratio of H2^CO/H2^CO throughout the 
disk. 
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Fig. 7. — The distribution of HCN and C2H2, over-plotted w ith isotherms (da shed lines) 
at the rotational temperatures of these molecules in GV Tau (IGibb et al.l 120071 ) - 115 and 
170 K, respectively. 
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f l2007h and lBrittain et al.l (120051 ) derived ^^CO/^CO ratios of 54±15 for GV Tau and 76±9 
for HL Tau. These ratios are typical of the transition region (GV Tau) and surface regions 
(HL Tau), assuming that these disks have an average ^^C/^^C ratio of 77. 



Gibb et al.l (120071 ) were able to calculate rotation temperatures for a number of species 



in GV Tau, including HCN (115±11K) and C2H2 (170±19 K), as well as ^^C O (200±40K) 
and ^^CO (260±20K). Using these and other constraints, iGibb et al.l (120071 ) were able to 
locate the origin of the hydrocarbon absorption to a region above the midplane, possibly 
the disk atmosphere, within ~10AU of the central star. These temperatures tally very 
well with our model: we have plotted the distributions of HCN and C2H2 from our model 
in Fig. [71 The c ontour lines show isotherms at the rotation temperatures calculated by 
Gibb et al.l (120071 ). The agreement is very good, especially given the very narrow vertical 
range in which these two molecules are abundant, due to freezeout in c older layers and either 
photo destruction or reaction with reactive ions in the layers above. IWhite fc Hillenbrand 
( I2OO4J ) make an estimate of the luminosity of GV Tau which is larger than that which we use 
in our model. Taking this into account might move the isotherm closer to the disk midplane, 
thus improving the agreement with the narrow molecule-rich layers. 



Gibb et al.l (120071 ) als o derived column den sities from their observations through the 



Brittain et all (l2005f ) calculated A^(i2CO)=7.5 x 10^^ cm'^ and 

and these column densities 



edge-on disk of GV Tau. 

A^(^^CO)=9.9 xlO^^cm"^ from their observations of HL Tau 



correspond to a radius of greater than 35 AU in our model. Derived rotational temperatures 
of 105 K and 80 K, respectively, indicate that the probed region is an intermediate layer of 
the disk. 



Methane abundances are also constrained by I Gibb et al.l (120071 . |200J), but rotational 
temperatures are not given. The gas temperature of the methane-rich layer in our disk is 
:^100K. In terms of column density, our calculated vertical column densities of iV(CH4)=3- 
21x 10^^ cm~^ at radii in the inner 10 AU of the disk ar e a factor of 20-300 above the observed 
upper limits for the GV Tau disk (IGibb et al.l 120071 ). This discrepancy could be down to 
differences in orientation between GV Tau and our model (column densities of methane 
and other molecules in our model are calculated vertically, whereas GV Tau is edge-on 
to the line of sight), or due to uncertainties in the chemistry of methane (perhaps in the 
acc urancy of activation e nergy barriers in the reactions which lead up to methane formation. 



see 



Agundez et al.l 120081 ). Abundances of methane are also dependent on initial conditions. 



with methane forming efficiently at the start of the interstellar cloud model on grains, due to 
high initial abundances of atomic C and H. Calculated column densities are given in Table 131 



Carr fc Najital (120081 ) have probed the inner few AU of the edge-on circumstellar disk of 



AA Tau, a typical classical T Tauri star, with observations from the Spitzer Space Telescope. 



- 29 - 



They report detections of the small molecules CO, C2H2, HCN, OH and H2O, with emission 
originating from within 2-3 AU of the star. The regions from which this molecular emission 
arises are hot (525-900 K). CO2 emission comes from cooler regions (350 K), and thus is 
assumed to come from larger radii. Our model shows very good correlation b etween regions 
of hig h abundance of these species and gas temperature with the results of ICarr fc Najita 
(120081 ). There is a slight disagreement with C2H2, which in our model is found in regions with 
tempe ratures of 70-250 K, significantly cooler than the 650 K calculated by 
( 2008 ). See Fig. [HI which shows distributions of the molecules detected by 



Carr fc Naiita 



Carr fc Najita 



(I2OO8I ) and isotherms at the designated temperatures attributed to those molecules. 



- 30 - 



Distribution of CO 




Radius [AU] 
Distribution of HCN 



1.5 - 




Radius [AU] 
Distribution of HjO 



1.5 - 




Radius [AU] 




12 3 4 5 

Radius [AU] 




12 3 4 5 

Radius [AU] 




Radius [AU] 



Fig. 8. — The distribution of CO, CO2, HCN, C2H2, H2O and OH over-plotted with isotherms 
(dash ed hnes) at the rotational temperatures of these molecules in AA Tau (jCarr fc Najita 
- 900, 350, 550, 500, 575 and 525 K, respectively. 
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Table 3. Selected vertical column densities of molecules in the model disk at 1, 10, 20 and 

30 AU 



opecies 


r^/^lm-riTi o^- 1 ATT 


Pnlnrnn of 1 D ATT 

v^oiumn dL lUAU 


Polnmn c»f 90 ATT 


Pnlnmn of 'in ATT 
V^OlUmn dL OU-ri-U 


H2 


8.3(24) 


1.2(24) 


5.7(23) 


3.7(23) 


He 


2.3(24) 


3.3(23) 


1.6(23) 


1.0(23) 


CO 


3.7(20) 


5.6(19) 


3.1 (19) 


2.4(19) 


CO2 


2.4(19) 


3.2(14) 


8.7(12) 


6.3(12) 


CH4 


2.1(18) 


2.5(17) 


1.8(17) 


4.3(18) 


C2H2 


1.2(19) 


3.2(15) 


7.0(13) 


4.9(13) 


C3H2 


6.0(18) 

V / 


1.4(16) 


8.7(15) 


7.7(14) 




1.6(18) 


3.3(12) 


4.3(11) 


2.9(12) 


C3H4 


6.6 19) 


2.5 16) 


9.2(11) 


5.0(12) 


C4H2 


2.1(18) 


8.1(14) 


1.1 (15) 


5.1 (14) 


H2O 


2.3(21) 


4.3(14) 


8.1 (13) 


7.2(12) 


N2 


1.6(19) 


2.1(17) 


1.1 (17) 


7.5(16) 


NH3 


2.3(20) 


1.2(15) 


3.7(13) 


2.0(12) 


HCN 


8.4(19) 


1.6(15) 


2.7(13) 


3.7(12) 


HC3N 


7.0(18) 


2.9(15) 


4.2(13) 


9.9(12) 




7.9(18) 


1.2(18) 


6.3(17) 


4.4(17) 


13CC2H4 


1.8(18) 


6.6 (14) 


2.8 (10) 


1.5(11) 


g-C2H2 


1.4(9) 


2.5(17) 


3.2(18) 


1.4(18) 


g-C2H6 


1.1(20) 


1.1 (19) 


8.4(17) 


2.9(17) 


g-C3H2 


2.2(8) 


3.1 (14) 


2.4(18) 


2.2(18) 


g-C3H4 


8.3(9) 


1.8(19) 


6.2(18) 


1.4(18) 


g-C4H2 


1.5(10) 


1.9(17) 


1.6(18) 


1.1(18) 


g-H20 


2.2(20) 


3.5 (20) 


1.7(20) 


1.0(20) 


g-NH3 


1.9(12) 


3.8(19) 


1.8(19) 


1.2(19) 


g-HCN 


1.2(15) 


9.7(18) 


4.6(18) 


3.0(18) 


g-'^CCHe 


1.9(18) 


2.0(17) 


1.8(16) 


6.1(15) 
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Note. — X (y) here represents x x 10^ cm ^. g- before a species denotes that that species is adsorbed 
onto a grain surface. 
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5.3. Significance for the Solar System 

Estimates of the mass of the protosolar disk suggest that it was somewhat more massive 
than the disk considered in this paper. To make a better comparison with Solar System data 
we have therefore also considered a model with M=10~^MQyr~^ and a=0.025, computed 
for us by Paola D'Alessio. This value of a results in a disk mass similar to the Minimum 



Mass Solar Nebula (MMSN; lHayashilll98ll ). All other parameters are the same as the fiducial 
model, apart from a new scaling for Zh, -2;h=0.0244R^^/^^. We input p^C]/[^''^C]=89 to the 
molecular cloud model to match present Solar System ^^C/^^C isotope ratios. 

Results from this model differ from the results of our fiducial model in two respects. 
Since the MMSN model is so much more dense, radial advection times are longer - a parcel 
of gas will pass from the outer edge of the model at 35 AU into the star in a time of 1.6 Myr, 
some four times greater than the fiducial model. This gives each parcel a greater processing 
time at each radius, which leads to: 

1) a greater range of degree of fractionation in the disk, such that, for instance, the 
fractionation in CO decreases from its input value into the disk of 50 to a minimum of 11, 
compared to 44 and 36, respectively, in the midplane of the fiducial model. 

2) a steeper fractionation gradient in the midplane. This increase in range of the frac- 
tionation occurs over the same distance in the model, thus the rate of change of fractionation 
is greater. This is also the case for vertical changes in fractionation, with often large changes 
in fractionation occurring over small distances (~1 AU). 

These differences can be seen in Fig. [Hlfor CO. 

Measurements of the ^^C/^^C ratio in the present day Solar System cluster around 
the telluric value of 89 (Fig. [TOj) . Although error bars in some cases can be very large, the 
majority of measurements are consistent with the idea that the bulk of Solar System material 
comes from a common origin, with an isotope ratio of 89. The most direct comparison 
we can draw between our MMSN model and present-day fractionation ratios is in the icy 
matter of comets, which is generally considered to be pristine. Cometary and meteoritic 
material is th ought to be remn ant from the very earliest phases of the formation of the Solar 



System (e.g., iMessengerl |2000| ). The degree to which this material has been processed is 
unknown, although clues can be found in analyses of interplanetary dust particles (IDPs) 
and carbonaceous chondrites (CCs) on meteorites. IDPs, which become trapped in the 
Earth's atmosphere, possibly originate in comets, whilst CCs come from the asteroid belt. 
Both types of compound may have been subject to heating, mixing and chemical reactions 
during the history of the solar system, possibly eradicating any chemical "history" . However, 
the results of carbon fractionation measurements show that in general, all measured comets 
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have a similar ratio (see Fig. [T0|) . Similarly, experiments replicating solar wind or cosmic-ray 
processing of comet surf aces has shown that these effects do not significantly contribute to 
carbon isotopic fraction (ILecluse et al.lll998l ). This is surprising, since ^^C/^^C ratios have 
been derived from both long- and short-period comets, which are thought to have formed in 
different regions of the protosolar disk. 

Carbon fractionation ratios in comets have thus far been determined by observing (one 
of) three molecules - HCN, CN and C2 - all members of the carbon isotope pool. CN is 
most likely a photodissociation product of HCN in cometary comae, although other parents 
may exist (lArpigny et al.ll2003l ). In our MMSN disk model CN and C2 ices are of very low 
abundance, which makes it difficult to give a ^^C/^^C ratio with a high degree of confidence. 
However, HCN ice is more plentiful, and has a ^^C/^^C ratio of 126-129 in the midplane 
region of the disk model. This ratio is inherited from the molecular cloud. It matches very 
well with observations of Comet Hale Bopp, a long-period comet: H^^CN/H^^CN = 110±12 



(iJewitt et al. 



1997 



Ziurvs et al.lll999h . but not with Comet Hyakutake: H^^CN/H^^qn 



34±12 (ILis et al.l 119971 ). To our knowledge these are the only three determinations of the 
H^^CN/H^^CN in comets. In the disk model, CN generally has a similar fractionation ratio to 
HCN. However, fractionation ratios derived from observations of CN in comets are somewhat 
different to HCN, in the region 65-115, with an average of ~90. This may be another 
indication of the alternative parentage of the CN molecule, or perhaps that there are some 
photo-fractionation effects in the photodissociation of HCN. 

What conclusions, then, can be drawn from this work given that our standard model of 
fractionation in a protoplanetary disk produces a very "stratified" result, where fractionation 
differs from region to region in the disk, and also from species to species? One possibility 
is that the protosolar disk was heated to some sufficiently high temperature for chemical 
exchange reactions to be ineffectual. This may have had the effect of "resetting" the carbon 
isotope ratios, similar to the process which occurred later in the formation of individual 
planets and the Sun. However, cooling from this state would have to have been relatively 
fast, faster than the chemical timescales for fractionation. 

A slightly different possibility relies on mixing material up to the surface layers of the 
disk, heating and photodestroying it, leading to a reset of the the carbon isotope ratio. For 
the processed material to be mixed down to the planetary accretion zones at the midplane of 
the disk would require vertical mixing timescales to be fast in relation to both the chemical 
timescales and accretion timescales. Given the high densities in the MMSN model, this 
seems unlikely since collision times and chemical reaction times will be short. 

A third consideration i s that of nebular shocks, which may have t ransiently heated ma- 
terial in the inner nebula (jPesch fc ConnoUyl l2002l : iKress et al.ll2002l ). Processed material 
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Fig. 9. — Carbon monoxide fractionation and distribution in the MMSN model. Compare 

to Fig. m 
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Fig. 10. — Measurements of the ^^C/^^C ratio in various objects of the Solar System. Filled 
circles indicate measurements of planets or the Sun and empty circles indicate measurements 
of planetary moons. Triangles indicate bulk isotope measurements of the ^^C/^^C ratio in 
meteorites, and have been placed at the radius of the asteroid belt. Comets, indicated by 
filled stars, have been placed outside of the radius of Neptune, for illustration, and similarly, 
IDPs (filled squares) have been placed at cometary radii to indicate their likely origin in 
comets. Please contact the author for a full list of references for the data in this plot. 
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could then have rapidly cooled and frozen out onto grain surfaces within chemical fraction- 
ation timescales, thus locking in the carbon isotope ratio of the hotter gas. However, the 
entire nebula would have had to have passed through such a shock to produce such a uniform 
carbon isotope ratio. 

So we conclude that whatever the mechanism which homogenised carbon isotope frac- 
tionation ratios in the Solar System (and these mechanisms bear future investigation), it 
occurred before comets and planets formed, yet after the initial collapse of the solar nebula's 
parent molecular cloud. 

6. Conclusions 

We have shown by means of a chemical model of a protoplanetary disk that the fraction- 
ation ratio of carbon, ^^C/^^C, varies according to position in the disk. The fractionation 
ratio is governed by temperature, which affects the rate and direction of chemical exchange 
reactions, and incident UV radiation, which affects self-shielding molecules such as CO. This 
produces a picture of the fractionation in a disk which is stratified. Certainly in the upper 
region of the disk, fractionation timescales are faster than mixing timescales, meaning that 
this stratified picture should persist if mixing were to be considered in the model. We also 
considered a Minimum Mass Solar Nebula model in which the chemistry in our infiowing gas 
packets has longer to evolve. Extremely good agreement was seen between observations of 
Solar System comets and the fractionation ratio in ices in the midplane of the disk model. 
In general, our chemical model has excellent agreement with recent observations of T Tauri 
disks, both in terms of chemical abundance and location in the disk. 
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A. Gas heating and cooling mechanisms 



A.l. Heating 

Photoelectric heating. Photoelectrons are emitted from dust grains when UV photons 
impinge upon the grain. The energy of these emitted particles depends on the energy of 
the incident photon and on the e mitting grain po t ential . The gas heating rate due to 
photoelectrons is approximat ed by iBakes fc Tielend (119941 ) , and is implemented as follows 
kamp fc van ZadelhofjboOlh . 



For the purposes of the thermal balance, we assume that there are two populations of 
dust grain within the disk - one consisting of silicate particles and the other consisting of 
a mixture of small graphitic and PAH grains. The silicate grains are assumed to have a 
radius, a = 0.1 /xm, and a UV c ross section per H nucleus (o"uv) of 5.856 x 10~^^ cm^/H 
atom (IKamp fc DuUemondl |200J) . The small graphite grains are assumed to have a size 
distribution ranging from 3-lOOAf lBakes fc TielensI Il994h . Lar ger grains do not contribute 
to the disk heating. 



Thus the photoelectric heating rate, FpE, is given by: 

TpE = bexiT-o ergcm"^s~^ 



(Al) 



( iBakes &: TielensI Il994j ) where e is an efficiency f actor dependin g on the material of the 
grain, x is the photon flux measured in units of the iHabing] (119681 ) fleld between 912-1110 A 
attenuated by dust and no is the total number density of hydrogen nuclei, ?T,o=2nH2+?T-H- 
6 is a constant and is equal to 1 x 10~^^ for graphite/PAH grains and 2.5 x 10~^o"uv for 
silicate grains, with auv = QabsTra^na, where Qabs is a UV absorption efficiency and na is 
the number density of dust grains. 



The photoelectric efficiencies and yields depend on material. Following iKamp &: van Zadelhofi 



(I2OOII ). we define a grain charge parameter, x = xT'^'^/ue, with He the electron density. Thus, 



6 X 10" 



epAH 



+ 



1.6 X lO-'ysiiT 



1.2 



1 + 1.8 X 10-%o-9i 1 + 1 X 10^2^ 
0.0487 5.8 X 10-^TO-^ 



+ 



1 + 4.0 X 10-%o-73 1 + 2 X 10-4 



X 



(A2) 
(A3) 



( IKamp fc van Zadelhofil200ll ). The yields for silicate material (ysn) are: 



Vsil 



0.70 for X < 10-4 
0.36 for 10-4 <x <1 
0.15 for a; > 1 



(A4) 
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Collisional de-excitation of H2. H2 molecules in excited rovibrational levels (due to 
Lyman- Werner band absorptions, for example) can decay to less energetic levels following 
collisions, thereby hea t ing th e gas in the process. This complex situation was simplified by 
Tielens &: HoUenbachl (119851 ) by considering a single excited pseudovibrational level of H2. 
They derived the following heating rate: 



CDx 



erg cm ^s ^, 



(A5) 



where nn is the number density of atomic hydrogen, tt-Hj is the number density of molecular 
hydrogen, tt-h* is the number density of vibrationally excited H2 (presumed to be a fixed 
proportion of H2, nHi=1 0~^^Ho) and E ^: is the effective energy of the pseudolevel (taken 



to be 4.166 x 10^^^ ergs, iLondoru Il978l ). 7^ and 7^^^ are the collisional de-excitation rate 
coefficients from the excited v=6 level to the v=0 level. Since these transitions tend to occur 
i n steps of Av=l, 7^ and 7^0 ' are assumed to equal to one sixth of the v=l-0 rate coefficients 
jTielens fc HoUenbachl Eossf l. Thus 7,^0 = 1.67 x 10-i3yTexp(-l 000 K/T) cm^ s'^ and 
2.33 X 10-^=^yTexp(-18 100 K/[T + 1 200 K]) cm^ s^^ 



7*0 



Photodissociation of H2. Whilst 90% of Lyman- Werner band excitation results in colli- 
sional de-excitation, the remaining 10% results in radiative diss ociation to a pair a hydroge n 
atoms, each carrying approximately 0.4 eV of kinetic energy (jStephens fc Dalgarndll973h . 
Hence the heating rate due to photodissociation of H2 can be approximated: 



Phd 



5.55 X 10 ^'^xrH2'^H2 erg cm s ^ 



(A6) 



( iTielens &: HoUenbachl Il985l ). Here F'jj^ is the photodissociation rate of H2, taking self- 
shielding into account (see Sect. I3.1."2l) . 



H2 formation. The formation and release of a hydrogen molecule from the surface of a 
dust grain contributes 4.48 eV of binding energy to the thermal balance. Since it is unclear 
how much o f that energy is transferred into rotat ion, vibration and translation , we follow the 
approach of iBlack &: Dalgarnd (119761 ) (and also iKamp fc van Zadelhofj I2OOII ) by assuming 
that one third goes into each motion. This results in a heating rate due to H2 formation: 



Form 



2.39 X 10~^^R 



form 



erg cm ^ s ^ . 



(A7) 



-Rform is the formation rate of H2 taken from lCazaux fc Tielend (j2002bl ). which incorporates 
revised H2 formation efficiencies (ena) at high temperatures, viz. 



form 



0.5n(H)t'H^idO"deH2'S'H cm ^s 



(A^ 
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a 



pc 



1 + 



(A9) 
(AlO) 



f H is the velocity of a hydrogen atom and Ud is the geometric cross section of a grain, is 
usually assumed to be ~10~^^?t,o, but can be calculated from the gas density: 



-1 



cm 



(All) 



We adopt a standard dust-to-gas mass ratio (5dg) of 0.01, a grain density (pd) of 2.5 gem ^ 
and a grain size (a) of 0.1/im. yU is the reduced mass of the gas (taken to be 2.4) and mn 



the mass of an H ato ni. Svt in Eg. lASlis the st icking coefficient of hydrog en and is assumed 
to be .4, in line with lCazaux fc Tielens The reader is referred to ICazaux Sz Tielend 

(|2002al E]) l or an explanation of the other terms. In the calculation of we assume F=10 
(Kamp, priv. comm.), where F is the accretion rate of H2 in units of monolayers per second. 
In the calculation of ^, we assume that when a grain is completely covered with a monolayer 
of ice no chemisorption of H atoms can occur, but physisorption can. Thus in Eq. lAlOt 
ii^Hc=-£'Hp=600 K, rather than i?Hc = 10 000K in the chemisorption case. 



C ionisation. Neutral carbon can be readily photoionise d in protopl anetary disks, with 



the release of approximately 1 eV of energy per ionisation (IBlack 



1987h. The heating rat e 



due to this process depends on the at tenuation by d ust absorption (IBlack &: Dalgarndll977l ). 
attenuation by self-a bsorption of C (IWernerlll970l ) and also attenuation by the H2 column 
( Ide Jong et al.lll980l ) (the first, second and third exponential terms, respectively). 



Clon 



2.2 X 10"^'n(C)xo 



exp(-2.4y4v - Tc - tu2V2/'JtvI) 
1 + ^H2^2/7rt;? 



erg cm 



(A12) 



(ITielens &: Hollenbachlll985l ). The various parameters are given by: 



TC = lO-i^iV(C), 
V2 = 9.2 X 10"^Sv"\ 



= 1.2 X 10-^^N(R2)Sv-\ 



(A13) 



Xo is the unattenuated photon fiux incident on the surface of the disk, A^(H2) and A^(C) are 
the column densities of molecular hydrogen and atomic carbon, respectively. 6v is the line 
broadening, taken to be equal to the sound speed. 
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Cosmic ray heating. Cos mic rays will penetrate the disk up until a certain column 



density of matter (150 gem , lUmebayashi &: Nakandll98ll ). and contribute to the thermal 



balance via the energy released in the formation an d subsequent electronic re combination of 
B.^. This process yields some 7eV per ionisation (iGlassgold &: Langeiill973l ). The heating 
rate is given by: 

Tc-ray = (1 + nHe/rio)CHno(1.28 X 10-" + 2.44 X 10-"(nHjno)) erg cm^^^ s"^ (A14) 

( IClavel et al.l Il978l ) , where nne is the number density of heliu m. For Ch, the cosm ic ray 
ionisation rate of hydrogen, we use the value of 5.98x 10"^^ s"^ (IWoodall et al.l 120071 ). 



X-ray heating. To calc ulate the heating effect due to X-rays we use the prescription of 
Gorti &: HoUenbachl (12004J ) for an X-ray energy spectrum of 0.5-10 keV: 

rx-ray= / exp(-iVoorx(^) )ax (^)no/hcatf^^ crgcm'^s'^ (A15) 



0.5 



HE) 



F{E) is the X-ray photon flux at a radius R: L(E) /AttR"^, where L(E ) is the stellar luminosity 
as a function of emitted energy (in keV). lGorti &: HoUenbachl (120041) fit a broken powe r law t o 
the X-ray spectrum of a weak-line T Tauri star presented in lFeigelson fc Montmerld (119991 ). 
resulting in: 

1.2Lx^-i-75 forE>2keV, 

0.18Lx^ forE<2keV. ^ ^ 

Lx is the X-ray luminosity of the central star in ergss^^, which is in general 10^'^-10-^L^ 
( jPeigelson fc Montmerldll999l ). We use Lx=10-'^ L^, and use the fit in Eq. IA161 since the X- 
ray spectrum of a classical T T auri star resembles that of a weak-line T Tauri star quite closely 
(IFeigelson fc Montmerld Il999l ). In Eq. IA15t A^o is the hydrogen nucleus column density of 
gas towards the central star, and crx{E) is the total X-ray photoabsorption cross-section per 
H nucleus: 



axiE) = 2.27 x 10-'"(E/1 keV) 



-2.485 



cm^/H nucleus. 



(A17) 



Jwilms et allboooh . /heat is the fraction of absorbed energy whi ch heats the gas, equal to 
0.1 for atomic gas and 0.4 for molecular gas (IMaloney et al.l 119961 ). 



There is also a secondary X-ray heating effect since the recombination of after 
ionisation results in a substantial population of H atoms in an excited state, which decay 
collisionally. Thus the primary rate (Eq. IA15P is augmented by an additional: 



r; 



X-ray 



2.22C|nH^2i ergcm-^s" 



(A18) 



( IShang et al.ll2002l ). is the X-ray ionisation rate of atomic hydrogen and E21 the energy 
gap between the n=2 and n=l excited levels of a hydrogen atom (10.19 eV). Cx is calculated 
according to the process described in Sect. 
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A. 2. Cooling 



The gas in a PPD is cooled by the line emission of atomic and molecular species, and 
also, when the gas temperature is greater than the dust temperature, by collisions between 
gas and dust. 

For the emission lines, we generally make use of the escape probability factor, /3(r). 
This factor takes into account the fact that atomic and molecular lines can have a much 
higher optical depth than that of the dust continuum. The escape probability factor, then, 
is the probability that a photon from a particular line with an optical depth, r, will escape 
from the disk. The maximum probability is 0.5, since we assume that any photons emitted 
in the negative z direction will be absorbed in the disk. 



Tielens fc HoUenbachl (119851 ) show that the optical depth in the vertical direction z is: 



8TTiy^{6v) Jo \nu[z')gi J 

where z=Q is the surface of the disk, contrary to our usual notation. In Eq. IA191 A^i is the 
Einstein transition probability from level u to level Z, v is the line frequency, 5v is the line 
broadening (assumed to be equal to the sound speed, cs), and rii are the level populations 
of levels u and /, and Qu and gi are the corresponding statistical weights. Ideally, the optical 
depth would be calculated by solving the level population equations in non-local thermal 
equilibrium (non-LTE). However, this is computationally expensive, and we have had to 
make assumptions in order to simplify the calculation. Hence we assume a plane-parallel 
slab, with uniform temperature and density. Level populations are in LTE. The optical 
depth then becomes: 



T[Z] 



exp 



kT 



- 1 



(A20) 



flTielens fc HoUenbachl 119851 ). 



The escape probability formalism is defined by Ide Jong et al.l (119801 ) : 

for r < 7, 

I 4. DOT 



1-cxp(-2.34t) 
4.68r 

N 0.5' 



4r In 



for r > 7. 



(A21) 



Atomic line cooling. To calculate electronic level populations, we solve the statistical 
equilibrium equations for a three-level system: 

2 

nj -Rjj = nj-Rjj and nx = ^^%, (A22) 
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where rix is the number density of species X and 



with 



AijP{Tij){l + Qij) + Cij for i > j, 
i9j/9i)Aji/3{^ji)Qji + Cji, for i < j, 



(A23) 



Qij Qji 



P{u,,) =B{uij, 2.7 K) + B{uij,Ta)Td{uij). 



(A24) 
(A25) 



where Cij is the colhsional rate from level i to level j and -P(^'jj) denotes the background 
radiation due to 1) the 2.7 K microwave background, and 2) the infrared emission of dust, ex- 
pressed in terms of Plan ck blackbody functions, with the dust temperature and r^(i/jj)=0.001 



(IHoUenbach et al.lll99ll ). Assuming that the local radiation field in the disk can be repre- 
sented by the sum of these two blackbodies is a simplification: it ignores the significant effect 
of stellar radiation in the disk surface and neglects the fact that in the disk midplane the 
local radiation is optically thick, and thus not proportional to Td(i'j.j ). This assumption will 
affect the cooling rate and thus the gas temperature (see, for example iKamp &: van Zadelhofi 
200l[ l. 



We consider the fine-structure lines of [Ol], [Cl] and [Cll], in the manner described above 
in Eqs. IA20HA24t calculating the statistical equilibrium in detail, since some regions of the 
disk are below the critical density required for LTE to be a reasonable assumption. 

The lowest three f ine structure lines for [Oil are at 63.2, 145.6 and 44.0 fim. Line data is 
taken from Table 1 of iKamp &: van Zadelhofa (l200ll ). who take into account collisions with 
H2, H and electrons. The cooling rate over all three lines is: 



Aoi = '^Pirui)huui[nu{0){Aui + BuiP{i^ui)) - ni{0)Bi^P{vui)] erg cm" 



■3s-i. 



(A26) 



(IKamp fc van Zadelhofa I2OOII ). h is the Planck constant, Vui the frequency of the fine- 
structure line, n„(0) and niiO) are the number densities of oxygen in the upper and lower 
levels of the transition, respectively, and Aui, Bui and Biu are the Einstein coefficients for 
the transition. 

A similar approach is taken for [Cl], although the critical density for the three [Cl] 
cooling lines at 609.2, 229 .9 and 369.0 /im is much lower, and LTE can be assumed. Again, 
atomic data is taken from lKamp fc van Zadelhofll (120011 ). The cooling rate is simply: 



A, 



CI 



/i[Aio/?(rio)z/io?^i(C)+A2o/?(T2o)«^2o^2(C) + A2i/3(r2i)z/2i?^2(C)] erg cm" 



(A27) 
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flKamp fc van ZadelhoflboOlf ). 

We also consider the [Cll] line at 157.7/im, again assuming LTE since the density in the 
disk is always above the critical density for this transition. The cooling rate is: 



A 



CII 



erg cm ^ s ^ , 



(A28) 



with atomic data from lKamp &: van Zadelhofj (l200ll ). 



Molecular line cooling. H2 rotational/vibrational and CO rotational line cooling are the 
main molecular cooling lines in the disk, although we have also taken into account cooling 
from CH molecules. 

The H2 cooling furi ction for the lowest 51 rovibrational energy levels has been derived by 
Le Bourlot et al.l (Il999l ). taking into account collisions with H, He and H2, but not including 
pumping by UV and X-ray radiation which will affect the H2 population in high energy 
levels. It is only strictly valid between temperatures of 10 ^-10'^ K and densities of 1-1 0^ cm~^. 
H2 cooling is negligible at densities less than lO^cm"^ ( iHoUenbach fc McKed 119791 ) and at 
temperatures less than 100 K, and we persist in using the function at the high densities which 
can be found in the disk. We adopt an ortho-to-para ratio of 1. 



26 rotational lines of CO are considered in calculating the cooling rate due to CO 



emission: 



25 



Ago = /5(rco) J] hiy^,[n,{CO){A.,, + B,,P{u,,)) - n,{CO)BjiP{ 



erg cm ^s ^. 



(A29) 



Collisional rate coeffici e nts a re from ISchinke et al.l (119851 ) . and other molecular data from 
Kamp &: van Zadelhofj (I2OOII ). The optical depth of CO lines is calculated using Eq. IA19I 



The emission lines of CH radicals have a minimal effect on the thermal balance in the 
inner disk, but ca n provide cooling in reg i ons w here CH is particularly abundant. We use 
the cooling rate of iKamp &: van Zadelhoffl (120011 ). 



AcH = (3{TcYL)noncnLr 



erg cm ^s ^, 



(A30) 



where the line coohng coefficient, Lrot, is derived by iHollenbach fc McKed (119791 ): 



^rot 



no-Eo(l+(ncr/no)+1.5Y^ ricr/no) 
kT{l-(n-ji^/no))atotV'j: 



erg cm s for Uq ^ rif. 



erg cm'^ s ^ for hq ^ 



(A31) 
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and /3(rcH)=l, which overestimates the coohng rate. The critical density for CH, ricr, is 
6.6xlOV(^/l 000 K) cm-3, and Aq and Eq are T.TxlO^^s"^ and 2.76x1 0-^^ erg respec- 
tively . The total inelastic cross-section of CH, atot, is taken to be Ix 10""^^ cm^ (IHoUenbach fc McKee 



19791 ). fx is the thermal velocity of the colliding hydrogen atoms. 



Cooling at high temperatures. At high temperatures (more than several hundred de- 
grees Kelvin) Lya and the line emission from the metastable ^D-^P transition of atomic 
oxy gen at 6300 A efScient l y coo l the gas. Cooling rates for both these processes are given 
by iSternberg &: Dalgarnd (1l989l ): 



Ao6300 



7.3 X 10"^^/?(rLy„)nenHexp(-118 400K/T) erg cm s 
1.8 X 10-2^(^06300)^0 exp(-22 800 K/T) x 



410n, 



(9T1/2 + 6.3 X lO-VT) 



erg cm s ^, 



(A32) 



(A33) 



with no the number density of oxygen. The third term in the brackets in Eq. IA33I takes into 
account excitation of the oxygen atoms by impacting electrons, as described in iDraine et al. 
(H). /?(rLya)=/5(ro63oo) = l. 



Gas-grain collisions. Collisions between dust grains and atoms and molecules will gener- 
ally act as an energy transfer mechanism, and will cool the gas if the gas temperature is larger 
than the dust temperature, which is generally the case in the regions under consideration. 
The heating rate due to these collisions is: 

Ag_g = 4.0 X 10"^2(7ra2)noradaT VT(r - Td) erg cm"^ s-\ (A34) 



where a=0.1 /im, the thermal accommodation coefficient, aT=0.3 (IBurke fc Hollenbachlll983l ) 
and nd is defined previously (Eq. lAlip . 
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